arXiv:1501.00205v2 [math.AP] 26 Jan 2015 


On the stability of some imaging fnnctionals 


Guillaume Bal * Olivier Pinaud ^ Lenya Ryzhik 1 
January 27, 2015 


Abstract 

This work is devoted to the stability/resolution analysis of several imaging functionals 
in complex environments. We consider both linear functionals in the wavefield as well as 
quadratic functionals based on wavefield correlations. Using simplified measurement set¬ 
tings and reduced functionals that retain the main features of functionals used in practice, 
we obtain optimal asymptotic estimates of the signal-to-noise ratios depending on the main 
physical parameters of the problem. We consider random media with possibly long-range 
dependence and with a correlation length that is less than or equal to the central wave¬ 
length of the source we aim to reconstruct. This corresponds to the wave propagation 
regimes of radiative transfer or homogenization. 


1 Introduction 

Imaging in complex media has a long history with many applications such as non-destructive 
testing, seismology, or underwater acoustics [iniESli. Standard methodologies consist in 
emitting a pulse in the heterogeneous medium and performing measurements of the wavefield 
or other relevant quantities at an array of detectors. Depending on the experimental setting, 
this may give access, for instance, to the backscattered wavefield (and its spatio-temporal 
correlations), or to the wave energy. The imaging procedure then typically amounts to back- 
propagating the measured data appropriately. When scattering of the wave by the medium is 
not too strong, measurements are usually migrated in a homogeneous medium neglecting the 
heterogeneities. This is the principle of the Kirchhoff migration and similar techniques, and is 
referred to as coherent imaging m- When scattering is stronger, the interaction between the 
wave and the medium cannot be neglected and a different model for the inversion is needed. 
We will only consider weak heterogeneities in this work, so that a homogeneous model is 
enough to obtain accurate reconstructions. We refer to mu for a consideration of stronger 
fluctuations and transport-based imaging. 

We are interested here in comparing two classes of imaging methods in terms of stability 
and resolution: one based on the wavefield measurements (and thus linear in the wavefields), 
such as Kirchhoff migration; and the other one based on correlations of the wavefield (and 
therefore quadratic in the wavefields), such as coherent interferometry imaging [13] • Here, 
“stability” refers to the stability of the reconstructions with respect to changes in the medium 
or to a measurement noise. 
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Stability and resolution of coherent imaging functionals in random waveguides were ad¬ 
dressed in m; for functionals based on topological derivatives, see also e.g. [2]. We refer 
to |19t 1201 [l2] and references therein for more details on wave propagation in random 

media and imaging. Our analysis is done in the framework of three-dimensional acoustic wave 
propagation in a random medium with correlation length ic and fluctuations amplitude of 
order cto- The quantities ic and cjo are fixed parameters of the experiments. The random 
medium is allowed to exhibit either short-range or long-range dependence. Our goal is to 
image a source centered at the origin from measurements performed over a detector array 
D. Denoting by A the central wavelength of the source, we will assume that it is larger than 
or equal to the correlation length, and that many wavelengths separate the source from the 
detector, so that we are in a high frequency regime. The case A = is referred to as the 
weak coupling regime in the literature (for an amplitude uo <C 1), while the case A 3> £c is 
the stochastic homogenization regime. The ratio £c/A is a critical parameter as it controls, 
along with do, the strength of the interaction of the sound waves with the heterogeneities. The 
regime A <C is addressed in [T] and is known as the random geometrical optics regime. The 
stability/resolution analysis is somewhat direct in that case, since a simple expression for the 
heterogeneous Green’s function is available by means of random travel times. This is not the 
case in our regimes of interest, where the interaction between the waves and the underlying 
medium is more difficult to describe mathematically. 

We will work with a simplified configuration and will define reduced imaging functionals 
that retain the main features of the functionals used in practice while offering more tractable 
computations. In such a setting, we obtain optimal estimates for the stability of the functional 
in terms of the most relevant physical parameters. We furthermore quantify the signal-to-noise 
ratio (SNR). The statistical stability is evaluated by computing the variance of the functionals 
at the source location. For the wave-based functional (WB), this involves the use of stationary 
phase techniques for oscillatory integrals, while computations related to the correlation-based 
functional (CB) involve averaging the scintillation function of Wigner transforms |23( 125] 
against singular test functions. It is now well-established that correlation-based functionals 
enjoy a better stability than wave-based functionals at the price of a lower resolution HH. 
This is due to self-averaging effects of Wigner transforms that we want to quantify here. 

The paper is organized as follows: in Sectionj^ we present the setting and our main results. 
We define the wave-based and correlation-based functionals and describe our models for the 
measurements in Section Section is devoted to the analysis of the resolution and stability 
of the correlation-based functional and Section [5] to that of the wave-based functional. The 
proofs of the main results are presented in Section and a conclusion is offered in Section 


2 Setting and main results 


2.1 Wave propagation and measurement setting 


The propagation of three-dimensional acoustic waves is described by the scalar wave equation 
for the pressure p: 

d‘^T) 

— = ■ [p~^{x)Vp], X E M^, t > 0, 

supplemented with initial conditions p{t = 0,x) and pt{t = 0,x). We suppose for simplicity 
that the density is constant and equal to one; p = 1- We also assume that the compressibility 
is random and takes the form 


k{x) = Kq 




Ko = 1 , 
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where (Tq measures the amplitude of the random fluctuations and £c their correlation length. 
We suppose that V is bounded and do is small enough so that k remains positive. The sound 
speed c = (k /?)“2 thus satisfies c = 1 + 0((Jq) ~ 1 and the average sound speed is cq = 1. The 
random field 1/ is a mean-zero stationary process with correlation function 

R(x)=E{V(x + y)V(y)}, 


where E denotes the ensemble average over realizations of the random medium. We assume R 
to be isotropic so that R(x) = i?(|x|). We will consider two types of correlation functions: 
(i) integrable R, which models random media with short-range correlations; and (ii) non- 
integrable R corresponding to media with long-range correlations. The latter media are of 
interest for instance when waves propagate through a turbulent atmosphere or the earth upper 
crust mm- Such properties can be translated into the power spectrum R, the Fourier 
transform of R, by supposing that R has the form 




= / e 


—ix’k 


R{x)dx, 


0<6 < 2 , 


where S' is a smooth function with fast decay at infinity. A simple dimensional analysis shows 
that R{x) behaves likes 5 < 3, as \x\ —)• oo, which is not integrable for large |x|. The 

case = 0 corresponds to integrable R since the function S is regular. We assume that 6 < 2 
so that the transport mean free path (see section 3.3) is well-defined. Propagation in media 
such that 2 < (5 < 3 is still possible but requires an elaborate theory of transport equations 
with highly singular scattering operators that are beyond the scope of this paper. 

Our setting of measurements is the following: we assume that measurements of the wave- 
field are performed on a three dimensional detector D and at a fixed time T. We assume that 
the detector Z) is a cube centered at xd = (T, 0,0) of side I < L; see figure 
then to image an initial condition 


2.1 Our goal is 


p{t = 0,x) = po{x), dtp{t = 0,x) = 0, 
localized around xq = (0, 0, 0). 

Practical settings of measurements usually involve recording in time of the wavefield on a 
surface detector. The two configurations share similar three-dimensional information: spatial 
3D for our setting, and spatial 2D -|- time ID for the standard configuration. Using wave 
propagation in a homogeneous medium, it is also relatively straightforward to pass from one 
measurement setting to the other. Our choice of a 3D detector was made because it offers 
slightly more tractable computations for the stability of the imaging functionals while quali¬ 
tatively preserving the same structure of data. 

As waves approximately propagate in a homogeneous medium with constant speed cq, it 
takes an average time to = {L — I j2)jco for the wave to reach the detector. We will therefore 
suppose that T > to, and even make the assumption that T = L/cq so that the wave packet 
has reached the center of the detector. Note that in such a measurement configuration, the 
range of the source is then known since the initial time of emission is given. We therefore 
focus only on the cross-range resolution. 

We consider an isotropic initial condition obtained by Fourier transform of a frequency 
profile g (a smooth function that decays rapidly at infinity): 

Po{x) = X^ [ g{B-Wk\ - \ko\/X))e-^^-^dk, 


where |A:o| is a non-dimensional parameter, A is the central wavelength, and Bcq ^ the band¬ 
width. 
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Rescaling variables as x — )■ xL, t —)■ tLjco, introducing the parameters rj = ^c/A, e = A/L, 
and still denoting by p the corresponding rescaled wavefield, the dimensionless wave equation 
now reads 

+ croR = Ap, p{t = (),x)=pl{x), ^(t = 0,x) = 0. (2.1) 

We quantify the bandwidth of the source in terms of the central frequency by setting B = 
{Lep)~^, where /r is a given non-dimensional parameter such that p S> 1. We suppose that 
p <C l/v/e so that the initial condition can be considered as broadband. More details about 
the latter condition will be given later. The rescaled initial condition then has the form 

pI{x) = g (^ep ^ (2.2) 

The normalization is chosen so that Pq(0) is of order 0(1). If p denotes the Fourier transform 
of g, we can write 

pUx) = Pl{\x\) c. IfcoP [ g ( — I 

Js2 y £p J 

Above, the symbol ~ means equality up to negligible terms and S'^ is the unit sphere of 
The initial condition is therefore essentially a function with support ep oscillating isotropically 
at a frequency {kol/e. As explained in the introduction, we assume that ic < A, which implies 


liX I 



Figure 2.1: Geometry 


g < 1. The case 77 ~ 1 leads to the radiative transfer regime in the limit e —)■ 0 127]. The case 
77 <C 1 corresponds to the homogenization regime [5| and a propagation in an effective medium 
(which is here the homogeneous medium since cq is small). The opposite case g ^ 1 gives rise 
to the random geometric regime addressed in [T]. 


For the asymptotic analysis, we recast the scalar wave equation (2.1) as a hrst-order hy¬ 
perbolic system on the wavefield u = (v,p), where v is the velocity: 


5v 

P^ + Vp = o, 


dp „ 


v = 0. 
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augmented with initial conditions p{t = 0,x) = Pq (x) and v(t = 0, x) = 0. The latter system 
is rewritten as 

(2,3) 


I + aoV 


(,i)) 


dt 


dxj 


where V = diag(0, V) G Dmn = 5m4<Jnj+<5n4<5mj is a 4x4 symmetric matrix for 1 < j < 3. 
Above and in the sequel, we use the Einstein convention of summation over repeated indices. 


2.2 Results 

Our main results concern the signal-to-noise ratio at the point x = 0 defined by 

SNR = I = or I^, 

aAMTKo) 

where stands for the CB functional and for the WB functional, which are defined 
further in section]^ We will also quantify the support of Var{/}(z). We will obtain optimal 
estimates for the CB and WB functionals in terms of the main parameters that define them 
as well as e, rj, 6 and p. Our measurements of the pressure are assumed to have the form 

p{t,x) = E{p}(t,x) + 6p^{t,x) + annp{t,x), 

where 6p^{t,x) models the statistical instabilities in the Born approximation and annp(t,x) 
is an additive mean zero noise with amplitude an- Note that by construction, coherent- 
based functionals will perform well only in the presence of relatively weak heterogeneities. 
Hence, modeling the statistical instabilities at first order is both practically relevant and 
mathematically feasible. Taking into account second order interactions is considerably more 
difficult mathematically; see [9] for a stability analysis of Wigner transforms (and not of the 
more complicated imaging functionals) in the paraxial approximation. The term cr„np models 
an additive noise at the detector and also takes into account in a very crude way the higher 
order interactions between the wave and the medium that are not included in 6p^. We suppose 
that Up oscillates at a frequency (so that a simple frequency analysis cannot separate the 
real signal from the noise) and that it is independent of the random medium. The variance can 
then be decomposed as Var{I^} = , where is the variance corresponding to the 

single scattering term and V)p the noise contribution. We also write Var{/^} = ■ 

Note that the measurement noise can actually be much larger than the average E{p}. Indeed, a 
simple analysis of E{p}(f, x) shows that it is of order e when |x| = 1 (omitting the absorption). 
It is the refocusing properties of the functionals that lead to a reconstructed source of order 
one from measurements of order e. 

We formally rescale the single scattering instabilities by e-coSt/2^ where is the transport 
mean free path defined in section so that the first-order interaction between the average 

field E{p} (proportional to and the medium has a comparable amplitude to E{p}. 

The fact that the single scattering instabilities are exponentially decreasing as the ballistic 
part can be proved in simplified regimes of propagation, where a closed-form equation for 
their variance can be obtained; see e.g. 017]. 

We need to introduce another important parameter, which is the typical length over which 
correlations are calculated in the CB functional. In dimensionless units, we define it in terms 
of the wavelength by Nqs. We assume that the detector cannot perform subwavelength mea¬ 
surements, which implies Nc > 1. We suppose for simplicity that correlations are calculated 
isotropically. Accounting for anisotropic correlations would add additional parameters and 
technicalities, but presents no conceptual difficulties. Since the resolution and the stability 


5 




of the CB functional are mostly influenced by Nqs and not by the size of the detector as 


by construction (see section 4.1 below), we will systematically suppose that the 
detector is sufficiently large compared to the wavelength so that its effects on the stability 
can be neglected in a first approximation. The parameter Nq is crucial in that it controls 


the resolution of the CB functional (shown in section 4.1 to be L/Nc) and its stability. Small 
values of Nc yield a good stability for a poor resolution, while large values lead to less stability 
with a resolution comparable to that of the WB functional. 


L 

I 

<70 

4 

(5 

A 

e, r] 
B 

<7n 

and 
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Nc 

AL 

I 

L 

Nc 


Table 2.1: N otati on 
Source-detector distance; see Fig. 2.1 

Size of the array; see Fig. |2.l| 

Amplitude of the random fluctuations 

Correlation length of the random fluctuations 

The correlation function decreases as 0 < <5 < 2 

Central wavelength of probing signal 

£ = X with L distance of propagation; f? = x 

Bcq^ is bandwidth with cq background sound speed 

/i = with width of probing signal equal to /rA 

Amplitude of detector noise annp{t,x) 

Correlation-based (CB) and wave-based (WB) imaging functionals, respectively 

Variance of the above imaging functionals 

eNc is length over which correlations are considered in 

Cross-range resolution of wave-based functional (Rayleigh formula) 

Cross-range resolution of correlation-based functional 


Inverse of transport mean free path; see (3.9) 


Our main results are the following: 

Theorem 2.1 Denote by V^{z) (resp. ) and V'^(z) (resp. ) the variances of the 
correlation-based and wave-based functionals defined in section\^for the single seattering con¬ 
tribution (resp. the noise contribution). Let be the transport mean free path defined in 


(3.9) in seetion 3.3 below. Then we have: 

3-5 


^-2SL 2 ,,-2 
V fUj ~ e h 


A 


V^(0) 


-2EL2 , 
~ e ^0 






KP(o) 

V^{0) ~ e 


CT„ 


L 


-C N, 


c 


-T.L2 
<70) 


V^i0)r^af. 


A {Ncfiy 

A (Ncur 


V 


^ N (3-5)A2 


N, 


c 
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Moreover, V^, and are mostly supported on \z\ < //e. Above, the notation ~ means 
equality up to multiplicative constants independent of the main parameters, and a A b = 
mm{a,b), a V b = max (a, 6) . Denote now by SNR^ (resp. SNR^) and SNR^ (resp. 
SNR^) the corresponding signal-to-noise ratios. Then: 

SNRgt (0) = SNR^ (0) A SNR^ (0) SNR^^ (0) = SNR^ (0) A SNR^ (0), 


where 


SNR^iO) ~ — 
(^0 


(3-5)/2 (((l-5)/2)A0)' 

"'-j (a) 


A 


({3-5)/2)A1 




SNR^iO)^ 

SNR^iO) 


e--2^^ 


(7n 




SNR^iO) 

<70 


-iSL 


e 2 


<7ri 


We below comment on the results of the theorem. 


Comparison CB-WB. Let us consider first the case of short-range correlations 5 = 0. We 
neglect absorption at this point (i.e. e~^^ = 1) and will take it into account when considering 
the SNR of the external noise. 

Let us start with the single scattering contributions SNR^ and SNR}^ . Suppose first 
that Nc is small enough and p large enough so that {L/NcX) A p = p. Then, 


SNR^ ~ — 
<70 



p. 


Hence, the SNR increases when the following quantities increase: \/ic (the dynamics gets 
closer to the homogenization regime where the homogenized solution is deterministic), p 
(smaller bandwidth, which results in a loss of range resolution) and (weaker fluctuations). 
The fact that the SNR increases as p stems from the self-averaging effects of the quadratic 
functional. There are no such effects for and we find SNR'^{0) = SNR^{0){pX/£c) so 
that the CB functional is more stable than the WB functional in this configuration of small Nc, 
whether in the radiative transfer regime (A = ic) or in the homogenization regime (A ^ ic)- 
It is shown in section 4.1 that the resolution of the CB functional is L/Nc, so that the 
best resolution is achieved for the largest possible value Nq, where correlations are calculated 
over the largest possible domain, namely Nc = l/X- Hence, the best resolution is XL/I, 
which is the celebrated Rayleigh formula and the same cross-range resolution as the WB 
functional. Now, for this best resolution and when Ip > L, then (L/NcX) A p = (L/NcX) so 
that SNR^{0) = SNR^ {0){X/ic) (equality here means equality up to multiplicative constants 
independent of the parameters A, Nc, £c, <7o)- This means that in the absence of external noise, 
and in a weak disorder regime where multiple scattering can be neglected, the SNR of the CB 
and WB functionals differ only by the factor {X/lc) for a similar resolution. This is a term 
of order one in the radiative transfer regime, and a large term in the homogenization regime. 
In the radiative transfer regime, in order to significantly increase the CB SNR compared 
to WB, one needs to decrease Nc and therefore to lower the resolution. This is the classical 
stability/resolution trade off as was already observed in pQ for the random geometrical regime. 
Note also that the statistical errors for both functionals are essentially localized on the support 
(of diameter pX) of the initial source. 
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We consider now the noise contributions and . A first important observation is that 
(x) is mostly localized on the support of the source, while (x) is essentially a constant 
everywhere. This stems from the fact that in the CB functional, the noise is correlated with 
the average field so that the functional keeps track of the source direction, while in the WB 
functional the noise is correlated with itself, see sections 3.3.2 and[^ Moreover,the SNR are 


SNR^ ( 0 ) 




(Jn 


SNR^(O) ~ SNR^{0) 




As for the single scattering SNR, for identical resolutions (i.e. when Nc = l/^), the SNR are 
comparable. In order to obtain to better SNR for CB, one needs to lower the resolution (and 
therefore decrease Nc) so that {L/{NcX)) Afi = ^ and SNR^{0) = SNR}^ (0)/r ;§> SNR^ (0). 


Minimal wavelength for a given SNR. We want to address here the question of the lowest 
central wavelength Am of the source (and therefore the best resolution) that can be achieved 
for a given SNR. In order to do so, we need to take into account the frequency-dependent 
absorption factor. We first consider the WB functional. We assume for concreteness that the 
noise contribution is larger than the single scattering one, but the same analysis holds for the 
reversed situation. Let us fix the SNR at one. Hence, 

Writing cr„ = e“ 
following expression 


vm has to be such that e = an 


this yields TiL = 2Tn- We will see in section 3.3.1 that S admits the 


e(l-l) A^-^(l-l)’ 


which gives 


A 


4-5 

m 


r\j 


7'm(l ~ |) 


Hence, as can be expected. Am decreases as the fluctuations of the random medium and the 
intensity of the noise decrease (i.e. do i and tO Remark as well that Am decreases as ^c 
does, and in a faster way compared to do or r„. This is also expected since the limit 0 
corresponds to the homogenization limit in which the wave propagates in a homogeneous 
medium provided do is small. The measurements are therefore primarily coherent and both 
the CB and WB functionals perform well. 

Let us consider now the CB functional and let <5 = 0. Since for {L/{NcX)) A /i ;§> 1, 
SNR^ is greater than SNR^, we may expect in principle to find a lower minimal wavelength. 
A way to exploit this fact is to consider a central wavelength of the source Am/a, for the 
Am above and some a > 1, and to compute correlations over a sufficiently small domain 
in order to gain stability, but a domain not too large so that the resolution is still better 
than Am- The resolution of the CB functional being L/Nc, we choose Nc = L/3/Am, with 
/? > 1 so that Xml 13 is smaller than Am- The prefactor in the SNR is then the square of 
{La/{NcXm)) A = {a/13) A ^ that we suppose is equal to (a//3)^. We find 


5iVR^(0) ~ d^“-i(a//3)2, 


and we look for a > 1 and (3 > 1 such that a'l^~^all3 = 1. Since d,i < 1, this is possible only 
when a is not too large (so that the central wavelength of the source cannot be too small, 
otherwise absorption is too strong) and when d„ is not too small either. When d^ is below the 
threshold, then the minimal wavelength is the same for the CB and WB functionals. Hence, 
compared to the WB functional, the averaging effects of the CB functional can be exploited in 
order to improve the optimal resolution only when the noise is significant. Finding the optimal 
Nc is actually a difficult problem, and is addressed numerically in |15) . 









Effects of long-range correlations. The strongest such effect is seen in S since S —)• oo 
as (5 —)■ 2. Hence, as the fluctuations get correlated over a larger and larger spatial range, the 
mean free path decreases and the amplitude of the signal becomes very small. This means 
that coherent-based functionals are therefore not efficient in random media with long-range 
correlations, and inversion methodologies based on transport equations mn that rely on 
uncoherent information should be used preferrably. 

There is also a loss of stability in the variance of the single scattering term for the CB 
functional, which for sufficiently long correlations (i.e. > 1) becomes larger than in the 

short-range case. In such a situation, for the SNR of the CB functional to be larger than 
that of the WB functional, one needs to decrease the resolution by a factor greater than in 
the short-range case. Notice moreover that there is an effect of the long-range dependence 
on the term that measures the distance to the homogenization regime. As the 

medium becomes correlated over larger distances, the variance increases, which suggests that 
the homogenization regime becomes less accurate. 


3 Imaging functionals and models 


We introduce in sections 3.1 and 3.2 the WB and CB functionals, and in section 3.3 our models 
for the measurements. 


3.1 Expression of the WB functional 

We give here the expression of the WB functional for generic solutions of the wave equation, 
which we denote by p and its time derivative dtp. The solution to the homogeneous wave 
equation with regular initial conditions {p{t = 0) = qo,dtp{t = 0) = qi), for {qo,qi) given, 
reads formally in three dimensions 

p{t) = dtG{t, ■)*qo + G{t, ■)*qi, G{t, x) = 7^<5o(i - kl), 

47r|x| 

where * denotes convolution in the spatial variables and 5o the Dirac measure at zero. In the 
Fourier space, this reduces to 

Tp{t, k) = cos \k\t J^qo{k) + Tqi{k), 

\k\ 

where we define the Fourier transform as Fp{k) = e~^^'^p{x)dx. 

From the data (p(T, x),dtp{T, x)) at a time t = T for x G D, the natural expression of the 
WB functional in our setting is obtained by backpropagating the measurements, similarly to 
the time reversal procedure [l8]. This leads to the functional 

I^{x) = [dtG{T, •) =. 1dp{T, •) - G(r, •) =. lDdtp{T, •)](x), 

where \d denotes the characteristic function of the detector. Using Fourier transforms, this 
can be recast as 

(^) = (cos \k\T H^dp{T, -m) - ^-^^^T{lDdtp{T,. 

This expression is slightly different from the classical Kirchhoff migration functional because 
of our different measurement setting. It nevertheless performs the same operation of back- 
propagation. Assuming that all the wavefield is measured, i.e. D = M^, and that = 0 so 
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Tp{t, k) = cos \k\t FqQ{k), one recovers the initial condition perfectly, i.e. (x) = qo{x). In 
practice, the entire waveheld is generally not available and diffraction effects limit the resolu¬ 
tion. In order to simplify the calculations, we assume without loss of generality that only the 
pressure p{T, x) is used for this functional. This modifies as 

I^{x) = (cos |A:|r T{1dp{T, •))(A:)) . (3.1) 

This signihcantly reduces the technicalities of our derivations while very little affecting the 
reconstructions. Indeed, for localized initial conditions, the value of the maximal peak of 
the functional is divided by a factor two compared to the full : if D = M^, and 
Tp{t, k) = cos \k\t Tqo{k), then 

((cos|A;|r)2.Tgo(A:)) = (-^9 o(A:)) + {^o^2\k\TTqo{k)) . 

The hrst term above yields ^qo{x) while the second one is essentially supported on a sphere 
of radius 2T far away from the source. We will analyse in section the expected value and 
the variance of the functional for random measured wavefields. 


3.2 Expression of the CB functional 


We define here the CB functional for a measured wavefield u. This requires us to introduce 
the Wigner transform of u [25l |23] and to decompose it into propagating and vortical modes. 
Since the initial velocity is identically zero, the amplitude of the latter modes remains zero at 
all times, see m- The projection onto the propagating mode is done with the eigenvectors 
of the dispersion matrix L{k) = kiD^ where D* was defined in (2.3). These vectors are given 
by b±{k) = [k, ±l)/\/2, with k = k/\k\ and we define in addition the matrices B± = b±'S> b±, 
where (8> denotes tensor product of vectors. The full matrix-valued Wigner transform of u is 
defined by 

W^{t, x,k) = / e" 


(27r) 




u(t,x-h y)dy. 


The quantity is a real-valued matrix. In our setting, the field u is measured at a time 
t = T = 1. We compute the correlations over a domain C C H, so that we do not form the 
full Wigner transform but only a smoothed version of it given by 


WI{T,x,k) = j^ [ e*"-^u(r,x-^)®u(r,x + ^)dy. 

(2vr)^ ee 2 2 

We only consider the propagating mode associated with the positive speed of propagation cq, 
the other one can be recovered by symmetry. This mode corresponds to the vector b+ and the 
associated amplitude is given by 


as{T,x,k) :=TT{Wl{T,x,k)fB+ik), 


(3.2) 


where {WfXy^ denotes the matrix transpose of IT| and Tr the matrix trace. The CB functional 
is then defined in our setting by 

{x) = / dk as{T, x + coTk, k), Dx = {k gM.^,{x + coTk) G D}. 

J Dx 


Above, we suppose implicitely that if (x-l-coTA:) G D, then C is such that (x + coTk) G D. 
The functional is slightly different from the classical coherent interferometry functional of [TT] 
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because of our measurement setting. However, the two functionals qualitatively perform the 
same operation, that is the backpropagation of field-correlations along rays. The difference 
resides in how these correlations are calculated: in our configuration, we have access to 3D 
volumic measurements at a fixed time so that the 3D spatial Wigner transform is available 
and is the retropropagated quantity; in [H], measurements are performed on a 2D surface 
and recorded in time so that the 3D spatial Wigner transform is replaced by a 2D spatio-lD 
temporal Wigner transform. 

For the stability and resolution analysis, it is convenient to recast in terms of the 
amplitude a associated with the full Wigner transform (i.e. computed for C = M^). We then 
obtain the expression, using the rescaled variables, cq = T = 1: 

I^{x) = [ f dkdqF^{k — q)a{l,x + k,q), F^{k) = f e^^'^dy. (3.3) 

JDx Jx^£yf2^G 

Note that is real-valued. 

3.3 Models for the measurements 

We introduce in this section our different models for the measurements. We will define a model 
for the mean of the measurements, and a model for the statistical instabilities. The latter is 
obtained by using the single scattering approximation (Born approximation). We start with 
the mean. 


3.3.1 Model for the mean 

The wavefield. Denoting by □ = 5^2 — A the d’Alembert operator and by 

psit) = dtG{t,-)*qo (3.4) 

the ballistic part associated to an initial condition go (with vanishing initial dtp), the solution 

d'^p 


to ( 2 . 1 ) reads: 


p{t,x) =PB{t,x) - CJoD 


-1 


V 


(t,x). 


^rjeJ dF 

Obtaining an expression for the expectation of p requires the analysis of the term EH((•/ {rje))d‘^ 2 P 
which is not straighforward and requires a diagrammatic expansion. Rather, we follow the 
simpler, heuristic, method of |24) . which amounts to iterating the above inversion procedure 
one more time and getting 


p{t,x) = 

PB{t,x) - croD"^ 


(3.5) 


H ( — 

rje 


d'^PB 

dF 


2rn-l 


(t, x) + 


92 


V { — \ ( — 


r]£ J dF 


d^p 
r]£ J dt'^ 


{t,x) 


and to replace p in the last term hy pB- Since E{p} ~ at first order in uq, pb may be 
replaced by E{p} to obtain a homogenized equation for E{p}. The result is [2Tj that a harmonic 
wave p{u},x) (Fourier transform in time of p) is damped exponentially, i.e. for |x| 7 ^ 0, 

1 — cos( 2 a;|x|) 


|E{:p(a;,x)}| < with 7 ( 0 ;) = 


-R I ^^ I dx. 


.7k 3 167r|a:|2 

The absorption coefficient 7 is obtained by adapting the setting of [23] to ours. Since our 
initial condition ( 2 . 2 ) is mostly localizec 
that the waves are absorbed by a factor 

l^o|^ _ o-qI^oP?? f 1 - cos( 27 |/co||x|) 


initial condition ( 2 . 2 ) is mostly localized around wavenumbers with frequency |fco|/e, we find 
7 1 ^ 1 = " , ' / - " "- RilxDdx := 7e. 


167r|x| 
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A Taylor expansion then leads to the classical \k\^ dependency of the absorption associated to 
the Rayleigh scattering: 


le ^ 


O-nlfcolV 

Stts 


R{\x\)dx = 


o-pIfeolV 

Svre 


R(0). 


We implicitely assumed above that -R(O) was defined, which holds in random media with short- 
range correlations but not in media with long-range correlations. The latter case is addressed 
below. We can relate the latter expression for 7 e to the mean free path ;= S“^(|A;o|) 
defined in the next paragraph in (|3.9[) by 


S = 


rf-al Trlfcol"^ 
e 2(27r)3 



R{r]{ko - \k\p))dp 




rj^al 7r|fco|^ 
e (27r)2 


R(0). 


Hence S ~ 27(|A:o|/e)- 

Obtaining such a relation in the case of long-range correlations is slightly more involved. 
To do so, we first notice that 


R 



(er/)3 *^5(0) 
(27r)3 



CiS{Q){£Vlf 


where the constant cs is given by C 5 = 2^ '^7r3/^r(^^)(r(|)) ^(27r) ^ [22], T being the Gamma 
function. The expression of 7 ^ is then 


7£~C5S’(0)7^ 


1 — cos 


2|fco||3;| 


167r|x 


5-5 


dx 


V ^r(^) I- cos(r) 

2'^7r3/2e 2r(|) Jo 


Recall that (5 < 2 so that the above integral is well-defined. The inverse of the mean free path 
is now in the long-range case 

^ 7^-Vg7r|fco|"-^5(0) f 1 ^ I 

£ 2(27r)32-5/2 Js2 (l-ko-py/^ ^ e (27r)32V2 (1 _ ^)5/2 

^ 5'(0) 

e 2^47 r(l-|)' 


We used again the fact that 5 < 2 to make sense of the integral. Since the following relation 
holds 

r(¥) r 1 - «s(r) , ^ 1 

0Fr(|)io 4(1 

we recover that S ~ 27 (|A:o|/e) in the long-range case. We will therefore systematically replace 
yd/col/e) by S/2 in the sequel. The expression of E{p}(t, x) is obtained by Fourier-transforming 
back E{p}(a;,x) and using the fact that |x| = cpt since the ballistic part is supported on the 
sphere of radius |x| = cot. This yields 

K{p}{t,x) ~ e~^^^^‘^pB{t,x). (3.6) 

This is our model for the average pressure. 
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The Wigner transform. It is shown in m that the Wigner transform of a wavefield 
u satisfies the following system 


dW^ 


Q 2 


V 2 


+ [Qi + —]W^ =\Vi + — ]W^ ■.= S^W^, 


with 


QiW 

ViW 


2 '^ 


00 

2 


dW dW . 

-1- 

dxj dxj 


V 2 W = i(TQ 


dydpe^P'y 

{2-kY 


dydpe^P'y 


V 


Q^W = ikjD^W - iWkjD^ 

x + ey\ ^-dWik + lp) 


er] 


dxn 


+ 


dW{k - \p) ■ / X + ey 


dxi 


D^V 


er] 


(2vr) 


V 


X + ey 
er] 


[k + p/2]jDm{k + p/2) 


-W{k - p/2) [k - p/2]jD^V 


X + ey 
er] 


Noticing that 


dy^p-yy 


X + ey 
er] 


= * E V{—pp), 


we can recast the r.h.s of (3.7) as 


S‘^W^ = Fe*pi 


e dW^ . £ 


+ 


edW^ . ,,, 


(3.7) 


" on *p F, 


Fe{x,p) = — (-) V(2r/p). 

e Vtt/ 

It is then well established m that a good approximation of E{Vh^} is 

E{WY-Wl:=Y,d%B±, 


where the matrices B± were introduced in section 3.2 and the amplitudes a/^{t,x,k) satisfy 
the following radiative transfer equation 




dt 


± cok ■ Vxa% = QiaY), 


(3.8) 


Q{a±){k) = / a{k,p){a%{p) - a%{k))5o{co\k\ - co\p\)dp, 

7r3 

n \ hr n 

a{k,p) = —^— ^^j^^R{p{k-p)). 


Above, (5o is the Dirac measure at zero. The equation (3.8) is supplemented with the initial 
condition 

a±(t = 0) = = Tr(W^’Y^±, 

where is the Wigner transform of the initial wavefield u(t = 0). Since al(A;) = a^(—/c), 
we focus only on the mode and drop both the + lower script and the e upper script for 
notational simplicity. We now need to compute the Wigner transform of the initial condition. 
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We consider an approximate expression that simplifies the analysis of the imaging functionals. 
The scalar Wigner transform of pg, denoted by Wq, reads: 


Woix,k) 


M 

( 2 vr) 


/52 




xg 


ki-{x- ey/2) \ - / ^2 • (x + gp/ 2 ) \ 




e/j 




Since the bandwidth parameter p is such that /je ^ e, we can separate the scales of the g terms 
and the oscillating exponentials. We can then state that, as is classical for Wigner transforms, 
different wavevectors lead in a weak sense to negligible contributions because of the highly 
oscillating term g^lf‘olx-(ki-k 2 )/£_ leading term in the initial condition is therefore 


WQ{x,k) ~ 


/52 




ki- {x- ey/2) \( h ■ {x + ey/2) \ 


eg, 




eg 




:= f dkiWi 

j 


ki 


X fc — I fco I fcl \ 


eg g 




where 

= ^3 dye'-^'^g {ki ■ {x - y/2)^ g (ki ■ {x + y/2)j . 

The effect of the function Wq in the phase space is essentially to localize the variable x around 
eg and the variable |A:| on a shell of radius |feo| with width g~^. This is what is expected 
since the initial condition for the wave equation is isotropic, oscillates at a frequency e~^\ko\ 
and has a bandwidth of order {eg)~^. The fact that g <C l/\/£ implies that the localization 
is stronger in the spatial variables than in the momentum variables, which can be seen as a 
broadband property. For the sake of simplicity, we replace the initial condition by a function 
that shares the same properties but has an easier expression to handle. We write: 


«o,± = gwo - l^ol)^ , 


where the function wq is smooth. With cq = 1, we then recast (3.8) as 

da 
dt 


+ k ■ Vxtt + S(/c)a = Q+(a), a±(t = 0,x, k) = gwQ —,^(1^1 - l^ol) 

eg 


(3.9) 


Q+{a){k) = [ a{k,p)a{p)5o{\k\ - \p\)dp, S(fe) = f R{g{k - \k\p))dp. 

JR3 e 2{27ry Js2 

When Uq = e, r/ = 1, we recover the usual equation for the weak coupling regime. Note that 
since R{x) = i?(|x|), we have R{k) = R{\k\) and S(fc) = S(|fc|). Going back to the measured 
amplitude defined in (3.2), an accurate approximation of its mean is therefore 

E{as} = TT{E{WI} fB+ = Tt{F- *k E{W^}fB+ ~ F" *k d, 


(3.10) 


where the filter was defined in (3.3). The integral solution to (3.9) reads, denoting by 
Ttd{x, k) := d{x — c^tk, k) the free transport semigroup: 

d{t,x,k) = TtaQ{x,k)e~^^^^'^^^ + f /c). (3.11) 

Jo 

The average d{t, x, k) is the sum of a ballistic term and the multiple scattering contribution. 
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3.3.2 Model for the random fluctuations in the measurements 


The wavefleld. We follow the Born approximation, which consists in retaining in (3.5) the 


terms at most linear in V. In order to take into account the absorption as explained in section 


2.2, we replace the ballistic term in (3.5) by IE{p}. In doing so, both the average IE{p} and the 


fluctuations are exponentially decreasing. The random instabilities on the pressure are then 
generated by the term 


6p‘^(t,x) =—(Toe 


-1 


V{ — 
r]£ 


d‘^PB 


{t,x) 


(3.12) 


We suppose that the additive noise on the wavefield u has the form cj„W(x) = cr,i(n^, n.p(x) = 
(T„(n^(|), np(|)) and is independent of the random medium. For simplicity, we assume that 
the noise entries are real and have the same correlation structure, E{n| (rE)n|(y)} = $(^), 
i,j = I,-- - ,4, where ‘h(x) = d>(|x|) and is smooth. Using (3.6) and (3.12), our model for the 
measurements is therefore 


p{t, X) = E{p}(t, x) + 5p^{t, x) + CTnUpix). 


(3.13) 


The Wigner transform. Let be the projection of the full Wigner transform onto the 
+ mode, i.e. a’^ = Tr(lF^)^B+, so that, for the as defined in (3.2), we have as = *ka^. We 

already know from (3.10) that Eja^} ~ *k a, with a the solution to the radiative transfer 
equation (3.9). We subsequently write a^ = a + 6a^, where 6a^ accounts for the random 
fluctuations. The simplest model for 6a^ is obtained for the single scattering approximation 
which consists in retaining in only terms at most linear in V, so that the related variance 
will be at most linear in R. This leads to defining the random term Wf 


dWl 


+ 


Qi + —^ Wf = 


e / 


0 ) 


(3.14) 


with vanishing initial conditions and where 


ot \ £ 


= 0 , 


w^{t = 0) = iu^’°. 


As was the case for the pressure, such a definition of Wf does not take into account the 
absorption factor We thus formally correct this and set 

Sa^ = e""°^*Tr(IUf)^5+. (3.15) 


The fluctuations of the amplitude 6a^ and the wavefleld can be related by writing the held u 
as E{u} + (5u, where Ju is obtained in the single scattering approximation. The perturbation 
5a^ is then the projection on the + mode of the sum of Wigner transforms lF[E{u},5u] + 
VF[(5u,E{u}]. Taking into account the external noise (T„n^(x) and denoting by af the projec¬ 
tion of VF[E{u},n‘'] + lU[n^,E{u}] on the -|- mode as in ( |3.15[ ), our complete model for the 
measurements in the single scattering approximation is therefore 


a^{t,x,k) = a{t,x,k) + 5a^{t,x,k) + anaf{t,x,k). 


(3.16) 


4 Analysis of the CB functional 

We compute in this section the mean and the variance of the CB functional. 
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4.1 Average of the functional 

Assume for the moment that our measurements have the form, for g a smooth function: 


as = *k g, 


g{x, k) = go{x - coTk, k) := Trgoix, k). 


Above, F^{k) is the hlter dehned in (3.3) and Tt is the free transport semigroup introduced 


in (3.11). Plugging this expression into the CB functional yields: 




'D^ 


dkdqF^{k - q)go{x + CoT{k - q),q). 


Let us verify first that when D = M^, we recover that I^[g]{x) = (27r)^ dqgo(x, q), so that 
if go is the scalar Wigner transform of a function ?/), we hnd I^[g\{x) = (27r)^|'0(x)p and the 
reconstruction is then perfect. This is immediate since from the definition of F^{k) given in 
(3.3) we conclude that F^{k) = (27r)^(5o(^) when C = When C is finite as is the case in 
practical situations, one does not recover dqg{x, q) but an approximate version of it limited 
by the resolution of the functional. This point is addressed below. 


Amplitude and resolution. For simplicity, we suppose that the domain C is chosen such 
that F^ defined in (3.3) is independent of x. We suppose in addition that C is a ball centered 


at zero of a certain radius. We parametrize C by 7 G [0,1] such that if G C, then we have 
|ey| < ro£^~'^ for some rg > 0. This means that the diameter of the ball is equal to 

roe“^ := Nc 


where the parameter Nc was introduced in section 2.2 Since the (rescaled) detector has a 


side ^ < 1, we have necessarily by construction that We could easily accommodate 

for anisotropic domains C with additional technicalities. In order to deal with a regular hlter, 
we smooth out the characteristic function of the unit disk and replace it by some approximate 
function x(x) = xd^^D- Rescaling y as ?/ —>■ yro£~'^, it comes for the hlter F^: 




/ 




ro\k\ 

£"< 


Recall that our measurements read as = F^ * o'", where is given by (3.16). The total 


functional is denoted by I^[a'^]{x) and its average is given by /‘^[a](x). The mean of the 
measurements d is the sum of a ballistic term a^{t,x,k) = TtaQ{x,k)e~^^^^^^^ and a multiple 
scattering term dehned in (3.11). Measurements in a homogeneous medium have the form 


TtaQ{x, k). The CB functional is tailored to backpropagate the data along the characteristics 
of the transport equation, so as to undo the effects of free transport. Since the multiple 
scattering term is smoother than the ballistic term, one can expect the inversion operation 
to produce a signal with a lower amplitude for the multiple scattering part than for the 
ballistic part. We therefore neglect the multiple scattering in the computation of the average 
functional. Moreover, the resolution of is mostly limited by the hlter F and not by the size 
of the detector since C is included in H. As a consequence, the limitations due to D on the 
resolution are negligible and we replace Dx by in the dehnition of . We then hnd. 


E{/C|a«]}(i) 




ge 


— 


£37 


'Da: 


ro\k - q\ 


£^ 


Wo 


X + k — q 

£g 


,K\q\ - l^ol) dkd\q\dq 
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so that the reconstructed image is essentially obtained by a convolution-type relation in the k 
variable between a localizing kernel and the source. This is clear when <C e/r where 


c. e-^/ [ F{\k\)wo dkd\q\dq. 

Jr3 Jks y J 

We can then define the resolution of the functional to be the scale of the localization, which 
is e^/ro = , or L/Nq in dimension variables. Hence, the larger the domain on which the 

correlations are calculated, the better is the resolution. When the domain is large, more phase 
information about the wavefield is retrieved, and this is the most useful information to achieve 
a good resolution. On the contrary, when 7 = 0 , the phase is lost and imaging is performed 
mostly using the singularity of the Green’s function, which lowers the resolution. 

Regarding the amplitude of the signal, we set x = 0 and compute the value of the peak. If 
the bandwidth parameter ^ with a > 1 — 7 , we find 

E{/^[a^]}(0) ~ f [ F {\k\)wo{0,\q\)dkdq e~^. 

JR3 Jr3 

In this case, the convolution is done at a smaller scale than the spatial support of wq, so 
that there is no geometrical loss of amplitude with respect to wo{x = 0). When a < 1 — 7 , 
convolution is done at a larger scale and we have 

E{/*^[a^]}(0) ~ / [ F {\k\) wq {eiOi + 6202,1^1) d\k\d9id92d\q\ 

jR'i Jr3 

Above, 9i and 02 are such that fe ~ g -t- Sfi{ei9i + 62 ^ 2 ) with ei ■ q = 62 ■ q = ei ■ 62 = 0 and 
|ei| = 1 62 1 = 1. The case 0 = 1 — 7 follows similarly. As a conclusion of this section, we 
therefore have 


4.2 Variance of the functional 


We compute in this section the variance of the CB functional with our measurements (3.16). 
For simplicity, we replace cq and T by their actual values, cq = T = 1. The total variance is 
defined by 

= E {I^[a^]\z)} - (E {I^[a%z)] f = E {l^[5a^]\z)] + a^E {l^[a%\\z)] 

:= V^{z) + V^{z). 

For some function ip to be defined later on, let us introduce the following quantity 
We{t) = / dxidx2dqidq2J{t,xi,qi,X2,q2)^{xi,qi)(p{x2.,q2)-, 

JK12 

where 

J = EjW^}, >V^(t,xi,g'i,X2,g'2) = 6a%t,xi,qi)5a^(t,X2,q2)- 

The function J is the single scattering approximation of the scintillation function of the Wigner 
transform and can be seen as a measure of the statistical instabilities. See for an exten¬ 

sive use of scintillation functions in the analysis of the statistical stability of wave propagation 
in random media. We have the following technical lemma, proved in section 6.1, that will be 
used to obtain an explicit expression for the variance V^: 
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Lemma 4.1 is given by 


Wsit) 




dkR{2'nk)\H^{t,k)\‘^, 


where the function reads 


Heit,k) = 


dsdxdqe^'^ £ (i^f + A:, s, x, q) 


/o 


-^1 'y ^ ^l^cr2fcFi,a2i ^2 k ■ q ^ ^ ^a2fcri,o 

(Ti,(T2=±l cri,(T2=±l 


„ - I A (X ~ <^isq 

fau<T 2 = h ■ Vx + t(T2\q\ I Wo ( - — -: 

(Pa 2 = ^{x + {t- s){q + a 2 k),q + a 2 k). 


hM - |A:o|) 


If we set 

^z{x,q) 

in the definition of then 



dkF^{k — q)So{x 


k — z) 


(4.2) 


V^{z) 


dpdqdkdk'F^{k — q)F^{k' — p)E{(Ia^(l, z + k, q)5a^{l, z + k',p)} 


dpdqdudvJ{l,u,q,v,p)pziu,q)pziv,p) = Weil). 


(4.3) 


The analysis of the variance is somewhat technical in that the scintillation is integrated 
against singular test functions (due to the application of the CB functional) while a classical 
stability analysis amounts to integrating against regular test functions. As before, we assume 
that the detector D is large enough so that its effects on the variance can be neglected and 
we replace it by as a first approximation. We have the following two propositions proved 
in sections 16.21 and 16.31 

Proposition 4.2 The variance of the CB functional for the single scattering contribution 
satisfies 

V^(z) ~ e- 2 ^c 72 r/ 2 A(i- 7 ) A V , |A < ep 

and is mostly supported on \z\ < pe. 

The above result gives the optimal dependency on the parameters p, s, 5, 7 and p. Regarding 
the contribution of the noise, we have: 

Proposition 4.3 The variance of the CB functional for the noise contribution satisfies 

Vniz) ~ for \z\ < ep, 

and Vff is mostly supported on \z\ < pe. 

As before, the result is optimal and the variance decreases as 7 goes to zero and resolution is 
lost. Notice that both variances are essentially supported on the support of the initial source. 
We now turn to the WB functional. 
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5 Analysis of the WB functional 

In this section we compute the mean and the variance of the WB functional. 


5.1 Average of the functional 


Using our model (3.13) and the definition of the functional (3.1), we have 




The cross-range resolution of is the same as the classical Kirchhoff functional and given 
by the Rayleigh formula XL/I. 


5.2 Variance of the functional 

The variance is defined by 

V^{z) = E{\I^\p]iz) - E {I^\p]} (z)|2} . 

As before, we distinguish the contribution of the single scattering, denoted by U^, from the 
one of the noise, denoted by . We then have the following propositions, proved in sections 
lOandlBAl 

Proposition 5.1 The variance of the WB functional for the single scattering contribution 
satisfies: 

{z) ~ e“^cJo for \z\ < e/r, 
and is mostly supported on \z\ < p,e. 

As for the CB functional, the results are optimal. Regarding the contribution of the noise, we 
have: 


Proposition 5.2 The variance of the WB functional for the noise contribution satisfies for 
all z: 





a 


2 

n' 


In the latter case, the variance is uniform in z and not just supported on the support of 
the source. This is due to the fact that, contrary to the CB functional which considers the 
correlation of the noise with the coherent wavefield (i.e. Vjp involves (p^)^(cj„np)^), the 
variance of the noise for the WB functional does not involve the average wavefield and the 
information about the source is lost (i.e. involves only (cnUp)'^). Remark also that 
and are comparable when 7 = 1, namely when both functionals have the same resolution. 
Vjp decreases when the resolution worsens. 


The proof of theorem 2.1 is then a straightforward consequence of (4.1), propositions 


4.3 5.1 


V/Arj. . J. lO Cli O Ui C4;J.^J.XUJ.Wi VV Cti. VJ. WJ. \-x.x 

and 5.2 after recasting e, g and rQe~'^ in terms of A, L, ic and Nc 
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6 Proofs 

6.1 Proof of lemma 14.11 


We start from (3.14) using the notations of section 3.3 We project the initial condition Wq 
onto the propagating modes and introduce the related amplitudes a± (recall there are no 
vortical modes since the initial velocity is zero): 

Wo(t,x,/c) = ^o±(t,x,/c)B±(/c), B± = h±®h± b± = ^(k,±lf', 


with 


da± 


± k • V,a± = 0, a±{t = 0) = = <(x/e, |A:| - |A:o|), 


where Wq(x, \k\ — |fco|) := ij wo{x/i j,, /i(|fc| — |fco|))- Split then in (3.7) into + 5| with 
obvious notation. We project (3.14) on the + mode and need to compute Tr(5fWg )^(A:)i?+(A:). 
Direct calculations yield 

Tv{SfWQ)'^{k)B+{k) = J dpfe{x, k-p){k-p±l) • Va;a± + i\p\a±^ 

/£(x,p) = — (-) e^V{2r]p). 
e Vtt/ 

In the same way, we find 

Ti{SIWqY' {k)B+{k) = ^ X] / ^Pfe{x, k-p){k-p±l) • V xa± - i\p\a^ 


and finally, recasting 5a^ by ai for simplicity of notation: 

9a 1 
Ik 


+ k ■ VxCLi = A^Se, p 1) • '^xa± + f|p|a±^ , 


Ae = ?fiJ dpfe{x,k-p)Se{x,p). 

Above, 3? stands for real part. The integral equation for ai reads: 

ai{t,x,k)= / dsAi;Ss{t — s,x — sk,k) := D~^Ai;Ss- 
Jo 

Let us introduce the product function 

W%t,xi,qi,X2,q2) ■= ai{t,xi,qi)ai{t,X2,q2)- 


We then have 

W^(t,Xi,gi,X2,g2) := {D~^A^Se){t,Xi,qi){D~^A^Se){t,X2,q2) 

_ 1 

“ 2 

X Se{t - Si,xi - siqi,pi)Seit - S2,X2 - S2q2, m)d'nid'n2dsids2 
1 

+ 2 

X Se{t - Si,xi - siqi,pi)Se{t - S2,X2 - S2q2,m)dr]idp2dsids2 

:=Ti+T2. 


3f? 



fe{xi - Siqi,qi - Pl)fe{x2 - S2q2,q2 - m) 


10 JO 


3? 



feixi - Siqi,qi - Pl)fe{x2 - S2q2,q2 - m) 


10 JO 
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Using the fact that E{U(^)U(z/)} = (27r)^i2(^)(5o(^ + ly), we find 


^{fe{x,p)fe{y,q)} = 
^{feix,p)feiy,q)} = 

Therefore 

X Seit - Si,Xi - Siqi,pi)Se 


R{2Tqp)5o{p + q) 

TT'^e^ 

aQlf 2ip-(:,-y) . 

-^e e R{2pp)do{p-q). 

77 -opz 


2i(,qi-Vl)-(xi-siqi-X2 + S2q2) 

r?i))e 

- S2,X2 - 5292, <?2 + 91“ r]i)dr]idsids2 


and 


T 2 


27r%2^, 



R{2p{qi - pi))e 


2i(qi-Vl)-(xi-siqi-X2+S2q2) 


X Seit - si,a:i - siqi,r]i)Seit - S2,X2 - 5292,92 - 9i + rii)dpidsids2. 


Hence 


E{W’^}(t,xi,g'i,X2,92) 


2e2 ^./o .0 


Ri2piqi - r]i))e 


2i(qi-Vl)-(xi-aiqi-X2 + S2q~2) 


Seit - Si,Xi - siqi,pi) 


X [Seit - S2,X2 - 5292, 92 + 9l “ "^l) + Seit - 52, X2 - 5292, 92 “ 9l + m)] 

2ik-{xi-sig\-X2+S2q2) 

III 

Jo Jo 


j Ri2pk)t 


X [Seit - si,xi - siqi,qi - k) + Seit - si,xi - siqi,qi + k)] 

X [Seit - S2,X2 - 5292, 92 + A:) + Seit - 52, ^2 - 5292,92 “ k)] . 


Accounting finally for the absorption e this allows to obtain the following expression for 

We: 

_2 3p-2S f 

= 332 / dkRi2pk)\Heit,k)\^ 

where 





2ik-x 

dsdxdqe ^ 


[5e(5,X,g 


/c) + Seis, X,q + A:)] pix + (t 


s)q,q). 


We have finally 


Seis,x,p) = J'^ik ■ p ±l)Jp ■ iV^w^) - |A:o|^ 

+ ^ ■ P ^ '^)\P\K \p\ - Ifcol^ , 

which replaced in He yields the expression given in the lemma. This is the final result and 
ends the proof. 
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6.2 Proof of proposition 4.2 


We compute formally the leading term in an asymptotic expansion of the variance, which will 
give us the optimal dependency in the main parameters. The proof involves the analysis of 
oscillating integrals coupled to localizing terms. We start from ( 4 . 3 ). We treat only the case 
7 G ( 0 , 1 ), the cases 7 = 0 and 7 = 1 follow by using the same approach. Plugging ( 4 . 2 ) in 
( 4 . 3 ), and assuming the detector is sufficiently large so that we can replace by in a first 
approximation, we find the following expression for H^-. 

HS = l,k) = Hl{k) + Hl{k) 


= 5 E 


o-i,cr2=ltl 


(Ti,(72=ltl 


dsdqduF.iu - q)k ■ 

4> 


0-2 


Sfl ' 




Efl 


where 


4 'g = 2; — ais{q — a2k) + sq — q + u 

fa2{x,q) = (^^q - cr2k ■ V X + i^J^cr2\q - cr2k\^ Wo{x, n{\q - a2k\ - |/co|))- 

The dependency of the function /0-2 on k and /r is not explicited in order to simplify already 
heavy notations. Expanding \H^ + leads to the following decomposition of w^: we can 
write We = wl + w'^ + wl '^, where 


^2 3„-2S 


dkR{2rik)\Hl{k)\‘^ 


and wl and wl'^ are obtained in the same fashion. We will focus only on the most technical 
term since lead either to the same leading order in terms of s, 7, rj and (Tq or to 

a negligible contribution. The term can itself be decomposed as 


= 


1 rl 


dsids2 (•••) + 


0 JO 


dsids2 (•••) + 


dsids2 (• • •) 


vl + V 2 + u|. 


The term can be shown to be negligible compared to the first two due to the mismatch of 
the integration domains, and we thus focus on vf and uf. We start with the most difficult 
term uf. We give a detailed derivation for the case fi = e“" with a < 1 — 7 (so that e/i <C 
and only state the result for the case a > 1 — 7. Let us start by denoting by ( 0 ^, 0 ^) and 
(0g,6lg) the angles defining u and q with the convention (0^,0^) £ (0, 27 r) x (0, vr), so that 
q = (cos sin6lg, sin6lq sin cos 0 g). After performing the changes of variables ( 0 u, 0 u) 
{0l,0g) + £fj{0u, 0u), we first find that 

u CE q + e^eiOl^F e^e20u with ei • g = 62 • g = ei • 62 = 0, |ei| = |e2| = 1. (6.1) 


Setting in addition si —>■ esi, S2 —)■ £82, \u\ —>■ |g| + e'>'|tt|/ro, k —>■ ek, .^1 = sigi and ^2 = S2q2, 
we find for the leading term, using the assumptions e/r <C <C 1 and /r ;§> 1: 


V2 


^2^3-5^4(l-7)+l-5^4^4g-2E 


/ / / 


dkd^id^2 




xG, 





( 6 . 2 ) 
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where 


GMli) = SiO) I dY\q\^Fi\u\i)f^^iz,, + (1 - ai)^ + e^Oi + e^Ol |g|, 1) 

f7i ,cr 2 =il 


dY = 




ei • ^ = 62 • ^ = ei • 62 = 0 


/° 2 (=^’| 9|>0 = y-^C-Yx + ifJ^cr 2 \q\j wo{x,n{\q\ - |A:o|)) 
and = z/efi. We will use the following relation: 


dk 


Ak-x 


\k\d 


{2tt)^6o{x) if /3 = 0, 
Ci 3 \x\^~^ if 0 < /? and 


/3/3,5,. 


(6.3) 


where Sq is the Dirac measure at zero and = 2^~ ^ 7 r 2 r(^-^)(r(^)) ^ (r being the gamma 
function). Note that (6.3) still holds for the case fd = 3,5, ••• provided the constant Cjs is 
adjusted, see [22]. We will not explicit this constant since only the dependency in x matters 


to us. Hence, we deduce from (6.3) that 


Ak-x 


dk^jj^kjk = 


5-1 


for 6 G [ 0 , 2 ], 


(6.4) 


where Cg is a constant. When 6 G (0,1), setting —>• s^i and ^2 —^ e ?2 in (6.2), it comes 

using (6.4) and /i ^ 1: 




^2^3-5^4(l-7)^4^4g-2E 


/ / 


d^id^2{^i ■ (6 - 6))(6 • (6 - 6 )) 

I 6 - 6 P-^I 6 PI 6 P 


G'.( 0 ,ei)G.( 0 , 6 ) 


plus a term of same order. The integral above is easily shown to be finite using the Hardy- 
Littlewood-Sobolev inequality [2^. A close look at Gz shows that only the term proportional 
to ^ ■ S/x is left because of the sum over a 2 - Setting finally |g| —>■ lA^ol + kl/l* yields ~ 

for 1^1 < fj,e. Because of the term it is not difficult to see that uf 


is mostly supported on |z| < fie. When 5 G [1,2) we use (6.4) and directly send e to zero in 


(6.2). The leading term in Gz is obtained for fii = 1, and we find 

f d^id^2 


^2^3-5^4(l-7)+l-5^4^4g-2E 


•^l£il<i •t| 6 l<i 1^1 
dY\q\^F{\u\Of^^{ze^ + ei9i + 62 ^^ \qli). 


G'°(|i)G 0 ( 6 ), 


where 

G°(O = S( 0 ) ^ 

CT 2 =dil 

Again, summation over a 2 in G^ imply that ~ for |z| < fie. The 

case 5 = 0 gives ~ and v^iz) is mostly supported on | 2 ;| < efi. Gathering 

all the previous results, we find for the case fi = e““ with a < 1 — 7 : 


vl{z) 


^2^3-5^4(l-7) + (l-5)A0^2^4g-2E^ 


|z| < efi, 


and V 2 (z) is essentially supported on | 2 :| < efi. When a > 1 — 7 (so that efi > e'^), very 
similar calculations show that for |z| < efi, we have Ws{z) ~ iTg 7 ^ '^ 6 ^^ ^'^'^^fi ^6 and that 
V 2 (z) <C t'KO) when \z\ S> efi. Combining the last two results, it finally follows that 

y2{z) ~ CTg773-'5e4(l-7)+(l-<5)A0^-2(^^^^) g7-l)4g-2S^ |^| < 
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and V 2 {z) is mostly supported on jzl < e^. 

We consider now the term uf. After the changes of variables si —)• esi and S 2 —)• SS 2 , we 
find 




dkR{2r]k)\V}{k)\^ 


where 


Vhk) 



(Tl,(72=ltl 


dsdqduF^{u 




z + u — q 
efi 



and / 0-2 is the same as before. Both cases a < 1 — 7 and a > 1 — 7 lead to the same order, and 
we therefore only detail the case a > 1 — 7 . The change of variables u —)• g + e'^ujro leads to: 


~ j j dsdqF{2e^ ^{{k-q)q-k)/{rQ\q\)))k-{q-a2k)e^^"^''ifa^{zef„q), 

(71,(T2=±1 


where F is the Fourier transform of F. Denoting then as before by {0\,9‘j^ and (0^,0^) the 
angles defining k and g, we perform the change of variables (0^, 0^) —)■ (0^, 0^) + roe^“'>'(0q, 0^), 
which yields 

g ~ fc + e^“'''roei0g + e^“"’'roe 20 g with ei • fc = 62 • ^ = ei • 62 = 0, |ei| = |e 2 | = 1. 
The term then becomes 

~ Y J J ds\Qfd\q\d9lde‘^ F{2\k\{ei9l + e29^)/\q\) 

(Tl,(72=ltl 

X sign(|g| - a 2 \k\)e^^^^^^ k\q\^ 


We need now to consider separately cj 2 = 1 and cj 2 = — 1 in the expression above. Let us 
start with <72 = 1, so that fa 2 =i involves the term wo{zs^, iJ,{\\q\ — |A:|| — |A:o|)), which suggests 
the change of variables |g| —)• \k\ + \ko\ + \q\/n- This brings a factor than cancels with 
the fi in the definition of / 0 - 2 . It remains to make sense of the integral in k in the definition 
of vf. Finite values pose no problems, and lead, after taking the square, to a term of overall 
order We focus therefore on large values of k. The integral of the term 

depending on s behaves as while the one in F is essentially independent of |A:| when k. 

Hence, taking the square, we find a term proportional to 


aoVo\3^"(i-7)g-2S 


'|fc|>i 


dkR{27jk)\k\ 


-2 


aoVo^3^"(l-7)g-2E 




d\k\S{2rik)\rik\ 


-5 


The integral is hnite when (5 > 1 and leads to a factor When 5 < 1, we set k —)• k/rj, which 
leads to a factor rf. The term associated to cj 2 = 1 is therefore of order 

Consider now the case cj 2 = —1, so that involves now wo{z£fj_, f^(lkl + |fc||-|^o|))- When 

\k\ > |/co|, the second argument in wq never vanishes, which leads to high order terms in 
The leading contribution is thus obtained for \k\ < |/co|- As mentioned above, bounded values 
of k leads to an overall order of which is negligible compared to the 

case 1 T 2 = 1 when 7 <C 1, or to the same order when r] = 1. It follows therefore that vf is of 
order and with the same arguments as uf, that it is mostly supported 

on l^l < e/r. Owing this, and taking the maximum with gives the result of the proposition. 
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6.3 Proof of proposition 4.3 


We use the same notation as in section 3.3.2 and only treat the case 7 G (0,1) for brevity. We 
have 

yn{z) = Ct2E{/^[<]2(z)} 

~ cj^ / dpdqdkdk'Fc{k — q)F£{k'— p)K{a'^{l, z + k, q)a^{l, z + k',p)}. 


Moreover 

aW^,x,q) = Tr(W[E{ri},n^] + W[n^E{n}])^(l,x,g)5+(g) 

= w[b+ ■ E{n}, b+ ■ n^] + w[b+ ■ n^, 5+ • E{u}] 

:= ^(u;[E{p},n^]+u;[n^,E{p}])+u;v, 

where tc[-, •] denotes the scalar Wigner transform and contains the remaining terms not 
included in the first one. We will not analyse the contributions to the variance of the terms 
involving Wv since they have the same structure as the one involving w\E{p},np ] + 'u;[n^,E{p}] 
and yield similar results. Then, 

^{ani'^,z,q)af,{l,z',p)} 

= k^{iw[lK{p},nl]+w[n‘p,E{p}]){z,q){w[E{p},nl]+w[n‘p,K{p}]){z',p)} + Re 


= ■^^{w[^{p},K]{z,p)w[E{p},ni]{z',q)} + R'^. 


(6.5) 


Again, we do not treat the term R'^ since it has the same structure as the first term above. 
We denote by V^p the contribution to the variance of the first term of (6.5). It reads, with 
the notation u = z + k and u' = z + k': 


T/C' ^2^-E 
^n,p ~ *^71® 


dpdqdkdk' dydy' Fe{k — q)Fe{k' — 

- ^v)pb{^,u - ~ ■ 


In Fourier variables and after setting and writing the cosines as sum of 

exponentials, this translates into: 


T/C' ^2^-S 

1/ rsj (T p 

^ n,p 


E 


0-l,(72=±l 


/ -u-e -u'e' 

/ dqdkdk'didi'Fe{k - q)Fe{k' + q-{^ + 

/e15 

(e - 2g) ,,(1^1 - \ko\)9pm - l^ol). 


Above, we used the fact that psit,^) = e^ 5 ^(£(|^| — |A:o|/e)) cost|.^| with g^ix) = pg{px). We 
then decompose ^ and as ^ = ^yA: + with ■ k = ■ k = 0. Denoting 

by and the angles defining k and k' and performing the changes of variables 

we find fc' ~ + eei0^, + £ 6201 , ei ■ k = 62 ■ k = 0 

and |ei| = \e 2 \ = 1. Together with ,^_i_ —>■ y/£^±, —>■ -v/eCj_; as well as |^y/c + \/eC±l — 
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1 ^ 11 1(1 + e|^_Lp/( 2 |^ll D), it comes, since <C e ^ and e <C e'^ as 7 G (0,1): 

V^^^r^ale^e-^ ^ J dXF,(A: - g)F,(fc|fc'| + g - ^(I^H | + l^f, |)/2) 

(Tl,(T2=ltl 

X g-*(ei®l,/+e26'^,)'(?||fc-2g)gi2-(^j-{^)/£ 

X $(^ 11 ^ - 2g)fi(^(|^||| - |fco|)5'/.(k||l - l^ol) 



Above, we used the notation = C\\k + y/£^± and + V^^±- The leading term is 

obtained for ^y +(Ti|^y| = +(T 2 |.^||| = 0 so that the first two phases vanish. There is otherwise 

some averaging that leads to a higher order contribution. Writing q = q\\k + q±, q± ■ k = 0, 
integrating in ( 0 ^, the phase term involving (ei 0 ^,+e 26 '^/)-(CyA:- 2 g) = - 2 ( 610 ^,+e 26 '^/)-Q', 
we obtain a Dirac measure enforcing that q± =0. As a consequence, using the fact that 
F{x) = T(|x|) and <I(A:) = ‘h(|fc|): 

~ E /T q\\\)Fe{\\k'\ ± gy - ^(|ey| + lefiDI) 

X e*l«ilV(2|«||l)g*l?llV(2|«f||)gi2.fc-€')/£ 

X ^(l?|| T2gy|)g^(|Cy| - |A:o|)5M(l?yI “ l^ol) 

[dY=[ [ [ [ [ [ [ dgyd^ydeyde±d^idM^'Pf- 

J Jm.+ Jr+ JR+ Jr3 Jr3 Jr+ 


The final expression is obtained by setting | A:| —)■ e'^| A;|/ro±(?y, (?y —)■ e'’'(?y/roT|fclT(|C|| l+l^y l)/2, 
l^yl —)• ^~^|^y| + |A:o| and |^||| — )■ + |^o|, which leads to, using that ^ S> 1: 







for some smooth function G and where Jq is the zero-th order Bessel function of the first kind. 
Hence, V^p(z) ~ iT^e^T- 7 ) 7 - 4 g-s fQj. and V^p{z) <C 1^(0) when \z\ » e/i so that 

V^p is essentially supported on \z\ < e/j.. This ends the proof. 


6.4 Proof of proposition 5.1| 

We use here the notation of section]^ With Sp^ the random fluctuations defined in (3.12), the 
variance of the WB functional admits the expression 


V^(z) = —E f dkdk'ed^-’^'k^ cos\k\cos\k'\E{(FlDSpn(t = l,k)(FlDSp^)it = l,k')}. 
{27ry Jjj3 

Similarly to the CB functional, we neglect the finite size of the detector as a first approximation, 
so that reads 


V^{z) ~ [ dkdk' cos |A:| cos \k'\E{{F6p^){t = 1 , k){F5jf){t = 1 , A;')}. 
JR3 
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We start by computing {F5p^){t,k). For a regular function /, we have 


ds 


(□ /) = / - s,-) * f{s,-)ds, so (J'(n f){t,k)= --sm\k\{t-s)I'f{s,k)ds. 

Jo Jo fI 


Using (3.4) and (3.12), this implies that 


{F6p^){t,k) = -doe |^sin|A;|(t-s)J'^U^—^ (s,A:) 

ft j„ ^ _ 

= -aoipefe~^^G — sm\k\{t - s){V{erj-) *k ^P^){s,k). 

Jo fI 

B _ 


Recall that p^ = dtG{t, •) *Pq, so that Ap^ = dtG{t, •) * Ap^ and 


l^ol 


Ap^{s, k) = -£'^\k\^g^, ( e(|A:|-) cos |A:|s, g^{x) = pg{px). 


Hence 


{F6p‘^){t, k) = aog^e^e 2 



0 4E3 


dsdp\k\ sin \k\{t —s)V{£pp)g^ ( £{\k — p\ — j cos \k—p\s 


Therefore, using that E{U(^)U(z 2 )} = {2tt)^R{^) 5o{^ + z^), we find for the second moment 


E{{F6p‘^){t = 1, k){F6p^){t = 1, k')} ~ 

(TQ(rye)^e® [ [ [ dsds'dp|A:||A;'| sin |/i;|(l — s) sin |fe'|(l — s')/(/c, s), 

Jo Jo Jm.3 

where 

f{k,k',p,s) = e~^R{epp)g^ (^{\k - p\ - 9jl - p\ - cos j/c-pls cos | A:'-p|s'. 

Hence, the variance of the functional reads 

V^{z) ~ f [ [ sin |/i;|(l — s) sin |/c'|(l — s') 

Jo Jo Jm.'^ 

cos |A;| cos |A:'|/(/c, fc',p, s). (6-6) 

Setting in order /c' —)• p + £~^k', k ^ p + £~^k, as well as p —>■ £~^p, and writing the sines and 
cosines as sum of complex exponentials, we find 


^ fJ3d4 J dX\k + p\\k'+p\h^{k',p,k) 

cl,a'2,0'3W4,0'5 W6=il 
2 Z 

X exp - {cJi|A; + p\ + (T2\k' +p|} exp - {a3\k + p\ + a^lk' + p\} 
z z 

X exp - {(dsl/sl - (T3IA: + p|)s}exp - {{aelk'l - a^lk' +p|)s'} , 

where 

h^{k',p,k) = e~^ed’^~’^'>^/^g^{\k\ - \ko\)gJl{\k'\ - \ko\)R{pp) 
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and 

dsds'dkdk'dp. 

The first two oscillating phases in compensate directly when cji = —as and a2 = —a^. 
The phases are otherwise strictly positive which leads to some averaging. The leading order 
is therefore obtained for ai = —as and a2 = —a4. We then use the following short time - long 
time decomposition of V^: 

V'^{z) = f [ dsds'{■■■)+ f [ dsds'{■■■)+ f f dsds {■■■)+ f f dsds'{---) 

Jo Jo Je Je Jo Je Je Jo 

■.= Y.V,iz). 

i=l 

The terms V3 and V4 can be shown to be negligible compared to the first two because of the 
different integration domains. We focus consequently only on the most interesting terms Vi 
since V2. Let us start with the most technical term V2. We perform a standard stationary phase 
analysis in the variable p. We are looking for a point po such that ask + pos + CJ4/C' + pos' = 0 
(first order term in the phase), together with as\k\ — as\k + po\ = ae\k'\ — a^lk' + po\ = 0 

(zero order term). This suggests that <75 = <73, ae = 174, po = 0 and asks = —a^k's'. Defining 

.^1 = sfc, ^2 = s'k', and performing the changes of variables p —)• ^/£p■, = —£73(74^2 + 

leads to, all computations done for the leading term: 

D2(-z) ~ asa^J W(2|fc'l))iaP 

<y3,iJ4=±l 

X k3/(2|A:|) +c74/(2|/c'|)|5|sin(6 -12)1^ 

X e-*("3-4|fc|+|fc'|)6-./e^^(|fc| _ |fco|)^(|fc'| - |A:o|) 

[ dX^ = [ [ f f d^id( 2 d\k\d\k'\. 

J Jr 3 Je<|6l<l J^+ J^+ 

Setting finally ^2 —^ £^2, \k\ —)■ |/co| + \k\/p, \k'\ —)■ |/co| + \k'\/p, we obtain 

V2iz) ~ aljf-^e-^ [ g 

Js^ 

The function above has a very similar structure to one of the initial condition, and is therefore 
mostly supported on \z\ < eg. 

Let us consider now the term Vi. After the changes of variables s —)■ ss and s' —)■ ss', we 
find (recall that we have for the leading order that ui = —CJ3 and £72 = —£74) 

Vi(z) ^ agp^ Y <^30-4 J dX\k + p\\k'+p\hi4{k',p,k) 

(T 3 ,(T 4 ,o- 5 ,( 76 =ltl 

X expf {(£75 |/c| - £73|A: +p|)s} expf {(£ 76 |A:'| - £74|A:' +p|)s'} . 

When 7 ] = 1 , then Vi is of the same order as V2- When p <C 1 , we will see that Li is the 
leading term. Remarking first that 

as\k+p\ expi{(£j5|A:| - £73!^ + p|)s} = £75!^! expi{(£75!A:| - as\k + p\)s} 

- -^expf {(f 75 |/c| - £ 73 |A: +p|)s} , 
i ds 
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with a similar observation for the other exponential, can be split into four terms. It can 
be shown that the leading one is the one involving the product of the derivatives, so that 


Vi(z) ^ f dkdk'dph^{k\p^k) 

cr3,(T4,cr5,cr6=±l 

X (expi{((T 5 |/c| -cr 3 \k + p\)} - 1) (exp i {(del A:'| - (J 4 ^\k' + p\)} - l) . 

With the changes of variables |fe| —>• |A:o| + \k\/p,, \k'\ —)■ |A:o| + \k'\/p and p —)■ p/rj, we find 


Vi{z) ~ doe' 




Js'2 Js^ 

V ) 

j dpR{p) ( 

' Q-io3\p\h) — ^ 


Terms in the expression above involving an oscillating exponential are negligible, which shows 
that Vi(0) ~ doe“^, and for the same reason as V 2 , Vi is mostly supported on \z\ < ep. The 
leading term is therefore Vi, which ends the proof. 


6.5 Proof of proposition 5.2 


The proof is straightforward, the variance of the WB functional for the noise contribution 
admits the expression 




cr„ 


_ '^n 

(27r)^ 


dkdk'e^^^ ^ cos |/c| cos \k'\E{{TlDnp){t = 1, k){Tl£)n^)it = 1, k')}. 


Neglecting the finite size of the detector in first approximation leads to 

{z) ~ [ dk{cos\k\f^{ek) ~ d^. 


7 Conclusion 

This work is concerned with the comparison in terms of resolution and stability of prototype 
wave-based and correlation-based imaging functionals. In the framework of 3D acoustic waves 
propagating in a random medium with possibly long-range correlations, we obtained optimal 
estimates of the variance and the SNR in terms of the main physical parameters of the problem. 
In the radiative transfer regime, we showed that for an identical cross-range resolution, the 
CB and WB functionals have a comparable SNR. The CB functional is shown to offer a better 
SNR provided the resolution is lowered, which is achieved by calculating correlations over a 
small domain compared to the detector. This is the classical stability/resolution trade-off. We 
obtained morever that the minimal central wavelength Am that the functionals could accurately 
reconstruct were identical in the regime of weak fluctuations in the random medium, and that 
in the case of larger fluctuations, the CB functional offered a better (smaller) Am (resolution) 
than the WB functional. 

We also investigated the effects of long-range correlations in the complex medium. We 
showed that coherent imaging became difficult to implement because the mean free path was 
very small and therefore the measured signal too weak compared to additive, external, noise. In 
such a case, transport-based imaging with lower resolution is a good alternative to wave-based 
(coherent) imaging. This will be investigated in more detail in future works. 
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